/*
 *  material.h
 *  EPPI-0.0
 *
 *  Created by Juan Carlos García Orden on 19/11/07.
 *
 *  Modified by Sergio Conde from March 2011 on.
 *
 */
#ifndef _material_h_
#define _material_h_

#include <iostream>
#include <fstream>
#include "commonutil.h"
#include <vector>
#include <math.h>
#include <MATH/TENSOR/tensor.h>
#include <MATH/TENSOR/xtensor.h>
#include <MATH/TENSOR/tutils.h>
#include <Eigen/Dense>

using namespace Eigen;

class Material
{
public:
	string name;
	string type;
	string hypertype;	 //Hyperelastic material type: Coupled or Decoupled ?

	// Scalar constitutive values

	double W;
	double Wiso;
	double U;
	double Ent;
	double Temp;

	// Stress tensor

	VectorXd S;
	VectorXd Svol;
	VectorXd Siso;
	VectorXd Qv;

	// Consistent Material tensor

	MatrixXd Ct;	// Linearized Constitutive tensor
	MatrixXd Ctvol; // Linearized Volumetric Constitutive tensor
	MatrixXd Ctiso; // Linearized Isochoric Constitutive tensor

    // Material properties
	double density;
    double E;
    double nu;
    double mu;
    double lambda;
    double K;

    // Invariants

	double I;	// First Invariant of the Right Cauchy-Green Tensor: I =  tr(C)
	double II;	// Second Invariant of the Right Cauchy-Green Tensor: II = 0.5 * (tr(C)^2 - tr(C^2))
	double III;	// Third Invariant of the Right Cauchy-Green Tensor: III = det(C)
	double J;	// Volume ratio J
	Vector3d Inv;

	double Im;	// First Invariant of the Modified Right Cauchy-Green Tensor: I =  tr(Cm), being Cm = J^(2/3)*C
	double IIm;	// Second Invariant of the Modified  Right Cauchy-Green Tensor: II = 0.5 * (tr(Cm)^2 - tr(Cm^2))
	double Jm;	// Modified volume ratio J^(-2/3);

	// Useful Tools

	Matrix3d CSym; 		// Symmetric Right Cauchy-Green Tensor in Matrix object
	Matrix3d CiSym;		// Inverse matrix of C
    Tensor2Sym Id, CTensor, CiTensor;
    Tensor4LSym Id4, IIc, P, Pt, tildeP;

public:
	Material();
	Material(const string&, const string &);
	~Material(){}

	// Hyperelastic material models

	virtual void calcW( const VectorXd & ){};

	virtual void calcS( const VectorXd & ){};
	virtual void calcSvol(const VectorXd & ){};
	virtual void calcSiso( const VectorXd & ){};

    virtual void calcCt(const VectorXd &){};

	// Thermoelastic/damage material models

	virtual void calcW( const VectorXd &, const double & ){};
	virtual void calcU( const VectorXd &, const double & ){};
	virtual void calcEntropy( const VectorXd &, const double & ){};
	virtual void calcTemp( const VectorXd &, const double & ){};
	virtual void calcConductivity( const double & ){};

    virtual void calcS( const VectorXd &, const double & ){};

    virtual void calcCt( const VectorXd &, const double & ){};

	// Viscoelastic material models
	virtual void calcW( const VectorXd &, const VectorXd & ){};

	virtual void calcS( const VectorXd &, const VectorXd & ){};
	virtual void calcQv( const int &, const VectorXd &, const VectorXd &){};

	virtual void calcCt( const VectorXd &, const VectorXd & ){};

	//Thermoviscoelastic material models
	virtual void calcW( const VectorXd &, const double &, const VectorXd & ){};
	virtual void calcU( const VectorXd &, const double &, const VectorXd & ){};
	virtual void calcEntropy( const VectorXd &, const double &, const VectorXd & ){};
	virtual void calcTemp( const VectorXd &, const double &, const VectorXd & ){};
    virtual void calcViscosity( const double & ){};

	virtual void calcS( const VectorXd &, const double &, const VectorXd & ){};
	virtual void calcQv( const int &, const VectorXd &, const double & , const VectorXd &){};

	virtual void calcCt( const VectorXd &, const double &, const VectorXd & ){};

	// Calculating useful tools
	virtual void calcInvariants( const VectorXd & ){};
	virtual void calcCTensor( const VectorXd & ){};
	virtual void calcCiTensor(const VectorXd & ){};
	virtual void calcP( const VectorXd & ){};
	virtual void calcPt( const VectorXd & ){};
	virtual void calctildeP( const VectorXd & ){};

	// Getting Strain Energy, Internal Energy, Entropy or Temperature
	double getW( void ){return W;};
    double getWiso( void ){return Wiso;};
    double getU( void ){return U;};
    double getEntropy( void ){return Ent;};
    double getTemp( void ){return Temp;};

	// Getting Stresses

	VectorXd getS( void ){return S;};
	VectorXd getSiso( void ){return Siso;};
	VectorXd getSvol( void ){return Svol;};
	VectorXd getQv( void ){return Qv;};

	// Getting Elastic tensor

	MatrixXd getCt( void ){return Ct;};
	MatrixXd getCtiso( void ){return Ctiso;};
	MatrixXd getCtvol( void ){return Ctvol;};

	// Second derivative of Strain Energy for 1D element (Trusses)
	virtual double getDDW( const VectorXd & ){};
	
    // Getting useful tools
	double getJ( void ){return J;};		// Volume ratio J. In a decoupled model it provides J^{-2/3}.
    Vector3d getInvariants( void ){Inv(0) = I; Inv(1) = II; Inv(2) = III; return Inv;}; // Provides the invariants of C
    Tensor2Sym getCTensor( void ){return CTensor;};  // Provides C in Tensor object
    Tensor2Sym getCiTensor( void ){return CiTensor;}; // Provides C^{-1} in Tensor object

    Tensor4LSym getP( void ){return P;};       // Provides the projection Tensor P
    Tensor4LSym getPt( void ){return Pt;};      // Provides the transpose of P
    Tensor4LSym gettildeP( void ){return tildeP;};  // Provides the modified projection Tensor

	// Getting material definition
	double getE(){return E;};
	double getNu(){return nu;};
	double getLambda(){return lambda;};
	double getMu(){return mu;};
	double getDensity(){return density;};
	double getBulkmodulus(){return K;};
	
	// Virtual functions for Composite
	virtual int getNfibers(){int n=1; return n;};

	// Virtual functions for (thermo)(visco)elastic case.
	virtual VectorXd getTau(){VectorXd t(2); return t;};
	virtual VectorXd getBeta(){VectorXd b(2); return b;};
	virtual VectorXd getEta(){VectorXd e(2); return e;};
    virtual VectorXd getVMu(){VectorXd m(2); return m;};
	virtual int getNIntVars( void ){int n=1; return n;};
	virtual bool getMode(){return true;};

    // Virtual functions for Thermoelastic
    virtual double getT0(){double T0=300; return T0;};          // Reference Temperature
    virtual double geta(){double a=1; return a;};               // Linear expansion coefficient
    virtual double getc0(){double c0 = 1.0; return c0;};        // Specific heat capacity
    virtual double getK0(){double K0 = 1.0; return K0;};        // Conductivity at T = T0
    virtual double getK(){double K = 1.0; return K;};           // Conductivity at T = T0
    virtual double getwk(){double wk = 1.0; return wk;};        // Softening parameter

    virtual VectorXd getViscosity(){VectorXd et(2); return et;};
};
#endif
